#### TABLE H.12: META-ANALYSIS NEGATIVE CONSEQUENCES
#### Table format of pooled effects for negative consequence outcomes

rm(list = ls())
source("./2_code/00_setup.R")

#### STUDY 1 ####

data <- fread(paste0(data_path, "data_study1.csv"), header = TRUE)
data$worked <- as.factor(data$worked)

# Welfare
lin_welfare_economic1 <- lm_lin(immigrant_welfare ~ treat_economic, covariates = ~
                                  edad + worked + male, data = data[(treat=='economic' | treat=='control')])
lin_welfare_humanitarian1 <- lm_lin(immigrant_welfare ~ treat_humanitarian, covariates = ~
                                      edad + worked + male, data = data[(treat=='humanitarian' | treat=='control')])

# Crime
lin_crime_economic1 <- lm_lin(immigrant_crime ~ treat_economic, covariates = ~
                                edad + worked + male, data = data[(treat=='economic' | treat=='control')])
lin_crime_humanitarian1 <- lm_lin(immigrant_crime ~ treat_humanitarian, covariates = ~
                                    edad + worked + male, data = data[(treat=='humanitarian' | treat=='control')])

# Employment
lin_employ_economic1 <- lm_lin(immigrant_employment ~ treat_economic, covariates = ~
                                 edad + worked + male, data = data[(treat=='economic' | treat=='control')])
lin_employ_humanitarian1 <- lm_lin(immigrant_employment ~ treat_humanitarian, covariates = ~
                                     edad + worked + male, data = data[(treat=='humanitarian' | treat=='control')])

# New Conflict Index
lin_conflict_economic1 <- lm_lin(new_conflict_s ~ treat_economic, covariates = ~
                                   edad + worked + male, data = data[(treat=='economic' | treat=='control')])
lin_conflict_humanitarian1 <- lm_lin(new_conflict_s ~ treat_humanitarian, covariates = ~
                                       edad + worked + male, data = data[(treat=='humanitarian' | treat=='control')])


#### STUDY 2 ####

data2 <- fread(paste0(data_path, "data_study2.csv"), header = TRUE)
data2$worked <- as.factor(data2$worked)

# Welfare
lin_welfare_economic2 <- lm_lin(immigrant_welfare ~ treat_economic, covariates = ~
                                  edad + worked + male, data = data2[(treat=='economic' | treat=='control')])
lin_welfare_humanitarian2 <- lm_lin(immigrant_welfare ~ treat_humanitarian, covariates = ~
                                      edad + worked + male, data = data2[(treat=='humanitarian' | treat=='control')])

# Crime
lin_crime_economic2 <- lm_lin(immigrant_crime ~ treat_economic, covariates = ~
                                edad + worked + male, data = data2[(treat=='economic' | treat=='control')])
lin_crime_humanitarian2 <- lm_lin(immigrant_crime ~ treat_humanitarian, covariates = ~
                                    edad + worked + male, data = data2[(treat=='humanitarian' | treat=='control')])

# Employment
lin_employ_economic2 <- lm_lin(immigrant_employment ~ treat_economic, covariates = ~
                                 edad + worked + male, data = data2[(treat=='economic' | treat=='control')])
lin_employ_humanitarian2 <- lm_lin(immigrant_employment ~ treat_humanitarian, covariates = ~
                                     edad + worked + male, data = data2[(treat=='humanitarian' | treat=='control')])

# New Conflict Index
lin_conflict_economic2 <- lm_lin(new_conflict_s ~ treat_economic, covariates = ~
                                   edad + worked + male, data = data2[(treat=='economic' | treat=='control')])
lin_conflict_humanitarian2 <- lm_lin(new_conflict_s ~ treat_humanitarian, covariates = ~
                                       edad + worked + male, data = data2[(treat=='humanitarian' | treat=='control')])


#### STUDY 3 ####

data3 <- fread(paste0(data_path, "data_study3.csv"), header = TRUE)
data3$employed <- as.factor(data3$employed)

# Welfare
lin_welfare_economic3 <- lm_lin(immigrant_welfare ~ treat_economic, covariates = ~
                                  edad + employed + male, data = data3[(treat=='economic' | treat=='control')])
lin_welfare_humanitarian3 <- lm_lin(immigrant_welfare ~ treat_humanitarian, covariates = ~
                                      edad + employed + male, data = data3[(treat=='humanitarian' | treat=='control')])

# Crime
lin_crime_economic3 <- lm_lin(immigrant_crime ~ treat_economic, covariates = ~
                                edad + employed + male, data = data3[(treat=='economic' | treat=='control')])
lin_crime_humanitarian3 <- lm_lin(immigrant_crime ~ treat_humanitarian, covariates = ~
                                    edad + employed + male, data = data3[(treat=='humanitarian' | treat=='control')])

# Employment
lin_employ_economic3 <- lm_lin(immigrant_employment ~ treat_economic, covariates = ~
                                 edad + employed + male, data = data3[(treat=='economic' | treat=='control')])
lin_employ_humanitarian3 <- lm_lin(immigrant_employment ~ treat_humanitarian, covariates = ~
                                     edad + employed + male, data = data3[(treat=='humanitarian' | treat=='control')])

# New Conflict Index
lin_conflict_economic3 <- lm_lin(new_conflict_s ~ treat_economic, covariates = ~
                                   edad + employed + male, data = data3[(treat=='economic' | treat=='control')])
lin_conflict_humanitarian3 <- lm_lin(new_conflict_s ~ treat_humanitarian, covariates = ~
                                       edad + employed + male, data = data3[(treat=='humanitarian' | treat=='control')])


#### CALCULATE OBSERVATIONS ####

# Economic treatment
n_obs_economic1 <- data %>%
  filter(treat == 'economic' | treat == 'control') %>%
  filter(complete.cases(immigrant_welfare, immigrant_employment, immigrant_crime,
                        new_conflict_s, treat_economic, edad, worked, male)) %>%
  summarize(count = n()) %>%
  pull(count)

n_obs_economic2 <- data2 %>%
  filter(treat == 'economic' | treat == 'control') %>%
  filter(complete.cases(immigrant_welfare, immigrant_employment, immigrant_crime,
                        new_conflict_s, treat_economic, edad, worked, male)) %>%
  summarize(count = n()) %>%
  pull(count)

n_obs_economic3 <- data3 %>%
  filter(treat == 'economic' | treat == 'control') %>%
  filter(complete.cases(immigrant_welfare, immigrant_employment, immigrant_crime,
                        new_conflict_s, treat_economic, edad, employed, male)) %>%
  summarize(count = n()) %>%
  pull(count)

obs_eco <- c(n_obs_economic1, n_obs_economic2, n_obs_economic3)

types_economic <- data.frame(
  treatment = rep('economic', 3),
  Study = c('Study 1', 'Study 2', 'Study 3'),
  n_total = obs_eco
)

# Humanitarian treatment
n_obs_human1 <- data %>%
  filter(treat == 'humanitarian' | treat == 'control') %>%
  filter(complete.cases(immigrant_welfare, immigrant_employment, immigrant_crime,
                        new_conflict_s, treat_humanitarian, edad, worked, male)) %>%
  summarize(count = n()) %>%
  pull(count)

n_obs_human2 <- data2 %>%
  filter(treat == 'humanitarian' | treat == 'control') %>%
  filter(complete.cases(immigrant_welfare, immigrant_employment, immigrant_crime,
                        new_conflict_s, treat_humanitarian, edad, worked, male)) %>%
  summarize(count = n()) %>%
  pull(count)

n_obs_human3 <- data3 %>%
  filter(treat == 'humanitarian' | treat == 'control') %>%
  filter(complete.cases(immigrant_welfare, immigrant_employment, immigrant_crime,
                        new_conflict_s, treat_humanitarian, edad, employed, male)) %>%
  summarize(count = n()) %>%
  pull(count)

obs_hum <- c(n_obs_human1, n_obs_human2, n_obs_human3)

types_humanitarian <- data.frame(
  treatment = rep('humanitarian', 3),
  Study = c('Study 1', 'Study 2', 'Study 3'),
  n_total = obs_hum
)


#### META-ANALYSIS ####

# Helper function to extract forest data
extract_forest_data <- function(meta_obj, narrative, outcome) {
  total_obs <- sum(meta_obj$n.e)
  data.frame(
    Study = c(meta_obj$studlab, "Pooled Effect"),
    Coefficient = c(meta_obj$TE, meta_obj$TE.random),
    Obs = c(meta_obj$n.e, total_obs),
    SE = c(meta_obj$seTE, meta_obj$seTE.random),
    Narrative = narrative,
    Outcome = outcome
  )
}


### Economic treatment, welfare outcome ###

out_economic_welfare <- rbind(
  tidy(lin_welfare_economic1) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_welfare_economic2) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_welfare_economic3) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error)
)
out_economic_welfare <- cbind(out_economic_welfare, types_economic)

m.gen_eco_welfare <- metagen(
  TE = estimate, seTE = std.error, studlab = Study,
  data = out_economic_welfare, common = FALSE, random = TRUE,
  method.tau = "REML", n.e = n_total, title = "economic_treatment_welfare"
)

forest_data_eco_welfare <- extract_forest_data(m.gen_eco_welfare, "Economic", "Welfare")


### Humanitarian treatment, welfare outcome ###

out_humanitarian_welfare <- rbind(
  tidy(lin_welfare_humanitarian1) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_welfare_humanitarian2) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_welfare_humanitarian3) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error)
)
out_humanitarian_welfare <- cbind(out_humanitarian_welfare, types_humanitarian)

m.gen_hum_welfare <- metagen(
  TE = estimate, seTE = std.error, studlab = Study,
  data = out_humanitarian_welfare, common = FALSE, random = TRUE,
  method.tau = "REML", n.e = n_total, title = "humanitarian_treatment_welfare"
)

forest_data_hum_welfare <- extract_forest_data(m.gen_hum_welfare, "Humanitarian", "Welfare")


### Economic treatment, crime outcome ###

out_economic_crime <- rbind(
  tidy(lin_crime_economic1) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_crime_economic2) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_crime_economic3) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error)
)
out_economic_crime <- cbind(out_economic_crime, types_economic)

m.gen_eco_crime <- metagen(
  TE = estimate, seTE = std.error, studlab = Study,
  data = out_economic_crime, common = FALSE, random = TRUE,
  method.tau = "REML", n.e = n_total, title = "economic_treatment_crime"
)

forest_data_eco_crime <- extract_forest_data(m.gen_eco_crime, "Economic", "Crime")


### Humanitarian treatment, crime outcome ###

out_humanitarian_crime <- rbind(
  tidy(lin_crime_humanitarian1) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_crime_humanitarian2) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_crime_humanitarian3) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error)
)
out_humanitarian_crime <- cbind(out_humanitarian_crime, types_humanitarian)

m.gen_hum_crime <- metagen(
  TE = estimate, seTE = std.error, studlab = Study,
  data = out_humanitarian_crime, common = FALSE, random = TRUE,
  method.tau = "REML", n.e = n_total, title = "humanitarian_treatment_crime"
)

forest_data_hum_crime <- extract_forest_data(m.gen_hum_crime, "Humanitarian", "Crime")


### Economic treatment, employment outcome ###

out_economic_employ <- rbind(
  tidy(lin_employ_economic1) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_employ_economic2) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_employ_economic3) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error)
)
out_economic_employ <- cbind(out_economic_employ, types_economic)

m.gen_eco_employ <- metagen(
  TE = estimate, seTE = std.error, studlab = Study,
  data = out_economic_employ, common = FALSE, random = TRUE,
  method.tau = "REML", n.e = n_total, title = "economic_treatment_employment"
)

forest_data_eco_employ <- extract_forest_data(m.gen_eco_employ, "Economic", "Employment")


### Humanitarian treatment, employment outcome ###

out_humanitarian_employ <- rbind(
  tidy(lin_employ_humanitarian1) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_employ_humanitarian2) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_employ_humanitarian3) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error)
)
out_humanitarian_employ <- cbind(out_humanitarian_employ, types_humanitarian)

m.gen_hum_employ <- metagen(
  TE = estimate, seTE = std.error, studlab = Study,
  data = out_humanitarian_employ, common = FALSE, random = TRUE,
  method.tau = "REML", n.e = n_total, title = "humanitarian_treatment_employment"
)

forest_data_hum_employ <- extract_forest_data(m.gen_hum_employ, "Humanitarian", "Employment")


### Economic treatment, new conflict index outcome ###

out_economic_conflict <- rbind(
  tidy(lin_conflict_economic1) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_conflict_economic2) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_conflict_economic3) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error)
)
out_economic_conflict <- cbind(out_economic_conflict, types_economic)

m.gen_eco_conflict <- metagen(
  TE = estimate, seTE = std.error, studlab = Study,
  data = out_economic_conflict, common = FALSE, random = TRUE,
  method.tau = "REML", n.e = n_total, title = "economic_treatment_conflict"
)

forest_data_eco_index <- extract_forest_data(m.gen_eco_conflict, "Economic", "Neg. Cons. Index")


### Humanitarian treatment, new conflict index outcome ###

out_humanitarian_conflict <- rbind(
  tidy(lin_conflict_humanitarian1) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_conflict_humanitarian2) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_conflict_humanitarian3) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error)
)
out_humanitarian_conflict <- cbind(out_humanitarian_conflict, types_humanitarian)

m.gen_hum_conflict <- metagen(
  TE = estimate, seTE = std.error, studlab = Study,
  data = out_humanitarian_conflict, common = FALSE, random = TRUE,
  method.tau = "REML", n.e = n_total, title = "humanitarian_treatment_conflict"
)

forest_data_hum_index <- extract_forest_data(m.gen_hum_conflict, "Humanitarian", "Neg. Cons. Index")


#### EXPORT TO LATEX ####

table_list <- list(
  forest_data_eco_welfare, forest_data_hum_welfare,
  forest_data_eco_crime, forest_data_hum_crime,
  forest_data_eco_employ, forest_data_hum_employ,
  forest_data_eco_index, forest_data_hum_index
)

dig <- c(0, 0, 3, 0, 3, 0, 0)

table_meta <- xtableList(
  table_list,
  caption = "Table Results of Meta Analysis",
  label = "table_metas",
  digits = dig
)

latex_output <- capture.output(print(table_meta, caption.placement = "top"))

writeLines(latex_output, paste0(tables_path, "table_H12.tex"))

